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We discuss a method of solving the time dependent Schrodinger equation for atoms with two active 
electrons in a strong laser field, which we used in a previous paper [ A. Scrinzi and B. Piraux, 
Phys. Rev. A 56, R13 (1997) ] to calculate ionization, double excitation and harmonic generation 
in Helium by short laser pulses. The method employs complex scaling and an expansion in an 
explicitly correlated basis. Convergence of the calculations is documented and error estimates are 
provided. The results for Helium at peak intensities up to 10 15 W/cm 2 and wave length 248 nm are 
accurate to at least 10 %. Similarly accurate calculations are presented for electron detachment and 
double excitation of the negative hydrogen ion. 

PACS numbers: 32.80.Rm, 32.80.Fb, 31.20.Di 
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I. INTRODUCTION 

Several programs are being pursued that aim at the 
description of three-dimensional two- or multi-electron 
atoms in strong laser fields The common moti- 

vation for these efforts is to obtain quantitative results 
for excitation, ionization, and generation of harmonics 
by laser pulses at intensities, where more than one elec- 
tron participates in the process. The various approaches 
emphasize different aspects of the problem. 

The first fully correlated three-dimensional calcula- 
tions for two-electron atoms in non-perturbative laser 
fields were made for constant laser intensity by the R- 
matrix Floquet method [ft] . Results have been published 
on H~ and He || and Mg ||. The advantage of the 
method is that it can be applied to multi-electron atoms, 
where existing atomic structure programs can be used. 
At large intensities, many angular momenta and Floquet 
blocks are required and very large systems of equations 
have to be solved. Recently an adaptation of the R- 
matrix method to solve the time-dependent Schrodinger 
equation has been proposed , which maintains the ap- 
plicability to general atoms, but may be less plagued by 
expansion size problems. 

The approach of Ref. ||] puts a strong emphasis on 
two-electron correlation in Tfe-like atoms at the expense 
of abandoning the realistic description of atomic struc- 
ture. The method solves the time dependent Schrodinger 
equation on a grid for the radial electron coordinates and 
with an expansion in single particle spherical harmonics 
for the angular degrees of freedom. By visualization of 
the wave function, in particular the process of direct dou- 
ble ionization could be studied. The implementation is 
adjusted to a massively parallel computer, but still the 
grid size and the length of the multipole expansion of the 
inter-electron potential is limited by computer resources. 

The method of Ref. Q describes a two-electron wave 
function by an expansion in numerical single electron 
wave functions that are calculated in a finite box. This 
provides a realistic representation of atomic structure and 
allows to adjust the wave function to the parameter range 



to be investigated. For example, when photoelectron 
spectra are to be extracted, continuous wave functions 
can be densely placed in range of electron energies of 
interest. Results have been published on Mg 

Here we present in detail the method employed for the 
calculation of excitation, ionization, and harmonic gen- 
eration in He published in a previous paper g]. Our 
purpose is to provide converged ab initio calculations for 
realistic laser parameters with special emphasis on elec- 
tron correlation. We use an expansion in explicitly cor- 
related two-electron basis functions and complex scaling 
H . The range of application is similar to that of Ref. [Q . 
The most important difference is our use of an explicitly 
correlated basis, which gives a very accurate description 
of atomic structure including doubly excited states with 
a relatively short expansion. The second crucial ingredi- 
ent of our method is complex scaling which, as we will 
show, gives a simple implementation of strictly outgoing 
wave boundary conditions by an L 2 method. The penalty 
of the method is the loss of a direct physical interpreta- 
tion of the continuous spectrum of the complex scaled 
operator. While this may not be a fundamental limita- 
tion of the complex scaling method, it does at present 
limit our results to total ionization, double excitation, 
and harmonic generation. 

Compared to Ref. f| we extend the calculations for 
Helium to higher laser intensities up to 10 15 W/cm 2 at the 
laser wave length of 248 nm. The foundation of the error 
estimates given in Ref. [|[ is presented and discussed in 
detail. By improvements of the basis the accuracy of the 
harmonic spectra calculation could be enhanced to about 
10%. We supplement the results by laser detachment and 
double excitation of H ~ . 



II. COMPUTATIONAL METHOD 

The Schrodinger equation of a two-electron atom ex- 
posed to a laser field described in velocity gauge with the 
dipole approximation is 
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with the atomic Hamiltonian 
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H = -^(Al + A 2 ) + j-r, — 

2 ri r 2 ri-r 2 
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where fi and r 2 denote the electron coordinates measured 
from the nucleus and Aj and are the corresponding 
Laplace and gradient operators. The nuclear charge is 
Z = 2 for Helium and Z = 1 for the negative hydrogen 
ion. Atomic units are used unless stated otherwise. 

The vector potential of a linearly polarized laser pulse 
is given by 

A{t) = h(t) sin(wi) (0, 0, A ) , (3) 

where we employed cos 2 and Gaussian shaped envelopes 

h cos 2(t) = [cos(7rf/T)] 2 
/iGauss(t) =exp[-(2i/T) 2 ] 

with the pulse duration T. 

The calculations were made in velocity gauge, since 
we found much better convergence than in length gauge, 
which is in agreement with previous experience and with 
theoretical arguments A calculation at intensity 

10 14 W / /cm 2 and frequencies of u> = 0.4 and 0.6 was re- 
peated in length gauge and gave the same results for 
ionization and single excitation. Results for double exci- 
tation could not be converged in length gauge. 

Eq. ([!]) is a 6 + 1-dimensional equation, which can be 
reduced to 5 + 1 dimensions because of cylindrical sym- 
metry, when the laser is linearly polarized. Due to the 
high dimensionality only a very limited range of the phase 
space can be numerically represented and one needs to 
carefully control the restrictions imposed on this space. 
The restrictions consist of basis set truncation and the 
boundary conditions at large distance. We first discuss 
the boundary conditions. 



A. Absorption of outgoing flux 

Ionization means that a finite portion of the wave func- 
tion moves away to arbitrarily large distances without 
further contributing to the dynamics of the system. In 
a finite space one needs to absorb this outgoing flux at 
the boundary of the space to avoid unphysical reflections. 
Common procedures are the use of a complex potential at 
large distances fll],|| or some form of mask function || . 
A more systematic control of the asymptotic boundary 
conditions has been proposed in Ref. [Q. Ideally, one 
admits only outgoing waves at large distances. However, 
outgoing wave boundary conditions are difficult to define 
in the presence of a dipole field, which ranges to arbi- 
trarily large distances. In any case, correctly imposed 



outgoing wave boundary conditions are energy depen- 
dent, which is, in general, quite difficult to implement 
computationally. An additional complication is that the 
resulting Hamiltonian is non-selfadjoint (the norm of the 
wave function on the finite space is not conserved) and 
it has non-orthogonal eigenfunctions. This may cause 
problems for computational implementations that rely 
on the orthogonality of the eigenfunctions of the Hamil- 
tonian. 

For calculations with only one active electron, which 
are effectively 2 + 1-dimensional, one can usually make 
the space large enough such that the boundary condi- 
tions are of secondary importance. In our case a more 
stringent method of absorbing outgoing flux is required. 
Such a method is complex scaling ||[l3]]. It consists in 
analytically continuing the Hamiltonian by multiplying 
the real coordinates by a complex number 



H{n ,r 2 ;t)^H 9 = H(e i0 n , e t6 f 2 ; t) , 



(4) 



where the scaling angle 9 is real and positive. For the 
time-independent Schrodinger equation, the mathemat- 
ical theory of complex scaling is well established. The 
new Hamiltonian Hg has the same bound state spec- 
trum as H, while the continuous spectrum is rotated 
by the angle —29 around the ionization thresholds into 
the lower half plane of complex energies. This separates 
the continua starting from different ionization thresh- 
olds. In the wedge-shaped area between the real axis 
and the rotated continua doubly excited states appear as 
square integrable eigenfunctions with complex eigenval- 
ues, whose imaginary parts give 1/2 of the autoioniza- 
tion widths. In an exact calculation, the values of bound 
state and resonance energies do not depend on the scal- 
ing angle 9. The method is being widely applied. For 
multi-photon physics it is used to calculate ionization 
rates and AC-Stark shifts of hydrogen-like systems by 
the Floquet method and in time-dependent calculations 
for hydrogen-like systems fl4|,|l5| ■ 

There is no complete mathematical theory for the ap- 
plication of complex scaling to time dependent problems, 
only partial results for the time evolution of bound and 
resonance states were found [[l6| . In the appendix we ar- 
gue that the restriction of the complex scaled Schrodinger 
equation 



8-$ 9 (fi,f 2 ;f) = Herein, r 2 ;t) 
at 



(5) 



to the space of square integrable functions is equivalent 
the unsealed equation with the constraint of strictly out- 
going wave boundary conditions. The outgoing wave so- 
lution at the coordinates (r*i,r 2 ) is obtained by evaluat- 
ing ^g at the back-scaled arguments (e~' ie ri, e _ie r 2 ). To 
establish this equivalence we need to assume far reaching 
analyticity properties of the solution r 2 ; t), which 

are difficult to prove in practice. 

Regardless of this mathematical problem, the method 
has been successfully employed in time-dependent calcu- 
lations pa and its validity could be verified numerically 
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. For hydrogen one can approach the limit 8 — > 0, i.e. 
directly compare with the usual Schrodinger equation. It 
was found that the projections on bound states and the 
expectation value of the dipole 

d(t) = (y e (e- w r,t)M*o(e- w r,t)) (6) 

do not depend on the scaling angle 8. The advantage of 
the complex scaled solution is that due to the absence of 
reflections a much shorter expansion of the wave function 
can be used when 8^0. For the two-electron system 
basis size requirements exclude very small scaling angles, 
but we found stable results for the excited state popula- 
tions and for d(t) in the range of 0.12 < 8 < 0.28 (see 
below) . 

The harmonic spectrum is obtained by Fourier trans- 
forming the acceleration of the dipole d. The total ion- 
ization yield is defined as 

iio„ = l-J2\(Mri, r 2 )\^e{e- i0 r u f 2 ;t = oo))| 2 , (7) 

i 

where $i is the i-th bound state function calculated with 
the real Hamiltonian Hg = o. We use the computationally 
more convenient formula 

lion = 1 - E I ^le (* >^)|^e(fi,f 2 ;t = oo))| 2 . (8) 

i 

where is the bound state eigenfunction of the com- 
plex scaled atomic Hamiltonian 

H ,e$i,e =Ei$ ifi . (9) 

Note the extra complex conjugation on the left hand 
function, i.e. in the integral the unconjugated function 
is used. Eqs. (]?]) and (||) are equivalent because of the 
analyticity of both, and *£g, and since for 8 = 
&i,e=o is real up to an overall phase (see appendix). 

The population of a doubly excited state a is deter- 
mined as 

P a = KK, d m\ 2 . (10) 

This equation does not have an unsealed analogue, since 
the resonance wave function <fr a ^g seizes to be square in- 
tegrable, when 8 approaches 0. 

B. Basis set expansion 

We approximate the solution of the complex scaled 
Schrodinger equation by expanding ^>g in a Hylleraas- 
like explicitly correlated basis 

L m „ L 
L=0 1=0 

k s m s n s 

X EE£ £ c w ;s (^ 2 m ^ 2 e- a » ri - /3sr2 (11) 

s fe=0m=0n=0 



The operator P\ projects on the singlet states and ri 2 := 
— The two-electron angular factors Gli for total 
angular momentum L and z-component L z = are 

Ga=r[rt l Y, C li,L-u- m yL{h)Y!:-\h). (12) 

m 

The Ci^. L _, _ are Clebsch-Gordan coefficients and 
are spherical harmonics. Note that for each L there are 
only L + 1 angular functions Gl.i- The major part of 
angular correlation, which in the usual atomic physics 
basis requires a large number of combinations of single- 
electron angular momenta I and L — I, is here contained 
in the inter-electron coordinate ri 2 . 

The expansion ( |ll| ) is known to be formally complete 
|Q and it converges rapidly for bound states of the 
three-body Coulomb system. In Ref. jl9) a further signif- 
icant improvement of the basis was achieved by selecting 
the combination of powers by the rule 

k + m + n + \k — m\(l — So n ) < p s . (13) 

This constraint can be understood as follows: The range 
of space covered by a basis function in the direction of 
r\ and ri is roughly k/a s and m/ 'f3 B , respectively. When 
a s ~ 0s and \k — m\ becomes large, the electrons remain 
far from each other and correlation, which is mostly con- 
tained in the coordinate ri 2 , becomes small. One there- 
fore needs fewer functions with ri 2 -dependence, when 
\k — m\ is large. The constraint leads to an important 
reduction in the expansion size without deteriorating the 
accuracy of bound and doubly excited state energies. 

In Refs. |ll],|2(]] for each state of Helium two sets of 
exponents were used, one describing the known asymp- 
totic behavior of the bound state wave function by se- 
lecting a\ — Z and /3i = V—Z 2 — 2E and a second 
one describing correlation by exponents a 2 = /3 2 , which 
were optimized to obtain the best bound state energy E. 
In our case we needed to describe many states, includ- 
ing strongly correlated doubly excited states, within the 
same basis set. Therefore we used several different sets 
of exponents for each L and I. As an example, the expo- 
nents and the powers used for L = 2 in the major part 
of the calculations are given in Table §. The first group 
of exponents is adjusted to describe the singly excited 
states and single electron continuum of the configura- 
tion type (Is, n'd), the middle group is for symmetrically 
doubly excited states and higher continua of the form 
(np, n'p) and the last group is for states (nd, n's) with 
single electron quantum numbers n, n' — 1,2,3. The 
particle-exchanged configurations are automatically in- 
cluded by the exchange symmetrization of the basis func- 
tions. 

It is important to observe that in velocity gauge double 
excitation must be included into the basis for a correct 
representation of the wave function, even when no real 
double excitation occurs. The reason is that the gauge 
transformation 
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*(fi,f 2 ;t) -» e^ (tH?1+ ^ ) *(f 1 ,f 2 ;i) 



(14) 



equally affects both coordinates and thus introduces vir- 
tual double excitation. 

Our basis functions Gli(ti, r 2 )r\r 2 n r\ 2 &~ OLari ^ sr ' 2 are 
strictly real, such that the phase as well as the 9- 
dependence of ^g(ri, r 2 ; t) is entirely contained in the 
expansion coefficients. This means that the expansion co- 
efficients of our complex scaled initial state $e,is(ri, r 2 ) 
are dependent on 9 and for each 9 and a different sys- 
tem of equations with different initial condition has to 
be solved. 

We can interpret the same fact in terms of the back- 
scaled solution ^e(e _te ri, e~ l V 2 ; t), which approximates 
the solution of the normal Schrodinger equation with 
outgoing-wave boundary conditions. The expansion 
functions for the outgoing-wave solution are then 

e- i ( L+k+m+n ^G u (r 1 ,r 2 )r%r%r? 2 e- e ~ iS ( a ' ri+0 ° r *\ 

i.e. they strongly depend on 9. By varying 9 we therefore 
vary the expansion functions for the physical solution 
and in this way we obtain an estimate of the basis set 
truncation error with respect to the radial coordinates 
rx, r 2 and r 12 . 



C. Alternative basis sets 



The good performance of expansion (11) for He is due 
to the fact that on the one hand it is very similar to 
the usual atomic physics expansion in terms of prod- 
ucts of single electron orbitals, which converges well for 
states, where the two electrons remain spatially sepa- 
rated. On the other hand the explicit dependence on ri 2 
allows a good description of the wave function at small 
inter-clcctronic distances, which is particularly impor- 
tant when both electrons are in the same shell, as in the 
ground state or in symmetrically excited states. 

We investigated several other expansions, which seem 
offer technical advantages or which are particularly suit- 
able for specific states of He. 

An implementation of arbitrary angular momentum 
for few-body systems is given by Wigner's _D-functions 
pl| , p2| . At high angular momenta that expansion allows 
a strong reduction of the number of non-zero matrix ele- 
ments in the calculation |23|. The D-functions separate 
the overall rotation of the system from internal degrees of 
freedom. The overall rotation is defined as the rotation 
between a body-fixed coordinate system, determined by 
two vectors a and b, and the laboratory coordinates. The 
D-functions carry indices D^ n , which designate the to- 
tal angular momentum L, the quantum number m of 
rotation around the lab-fixed z-axis, and the quantum 
number n of the rotation around a. Although there is 
a certain freedom of choice for a and b, they cannot be 
identified with the electron coordinates f\ and r 2 , since 



the definition of the D-functions is not symmetric un- 
der exchange of a and b. In order to implement electron 
exchange symmetry, one can for example use the Jacobi 
coordinates 



a = r\ 
b = r i 



r 2 
f 2 



as the body-fixed vectors. With this choice, the sub- 
scripts of D^ nn refer to collective rotations of the elec- 
trons. This is desirable for some highly correlated states, 
like the Wannier states . In the unsymmetrically ex- 
cited states that dominate the wave function of an atom 
excited by a laser pulse, one electron carries the major 
part of angular momentum (except for symmetrization) , 
which leads to poor convergence of the Z?-function ex- 
pansion. 

Similar problems arise for an expansion with respect 
to the perimetric coordinates 



u = — n - 
v = r*i 
w = r% 



r 2 + ri2 

- r 2 + r 12 

- r 2 - r 12 . 



While the interparticle coordinates r\,r 2 , r± 2 are subject 
to the triangular inequality \r\ — r 2 \ < r\ 2 < r\ + r 2 , 
the perimetric coordinates u, v, and w each vary inde- 
pendently in [0, oo). This simplifies the calculation of in- 
tegrals and allows to find expansion functions, where the 
operator matrices become sparse. Like with the Jacobi 
coordinates, the unsymmetrically excited states, where 
the two electrons move largely independently, are not ef- 
ficiently described by such an expansion, since u, v, and 
uu each contain both coordinates r\ and r 2 . It also ap- 
pears difficult to find a constraint like ([l3|) to cut down 
on the basis size. 

Finally, we explored an expansion with respect to the 
coordinates r±, r 2 , and cos#i 2 := fi-f 2 . This is very sim- 
ilar to an expansion in single electron orbitals. The main 
advantage is that the calculation of matrix elements be- 
comes simple. However, the expansion length is generally 
larger than with the explicitly correlated basis and bound 
state accuracies beyond 10~ 4 a.u. become extremely hard 
to achieve fE3|. 



D. Numerics and computation 

The expansion ( |ll| ) is notoriously numerically difficult. 
The main reason is that the metrical matrix 

Sa ■■= (15) 

for the basis functions 

I*) = G u (f u f 2 )r k 1 i r^e- a ^-' 3 ^ 



(16) 



rapidly becomes ill-conditioned with increasing powers 
ki, mi, rii. The situation is further aggravated by our 



4 



use of several different sets of exponents a s , f3 s for the 
same L, which makes the basis formally overcomplete. 
We were able to control these problems by performing 
accurate integrations, by appropriately normalizing the 
basis, and by removing near singular values from the met- 
rical matrix. Still, at high angular momenta we needed 
to resort to Fortran REAL* 16 (« 32 decimal digits) ac- 
curacy in the calculation of the matrix elements. 
We first express the angular factors in the form 



r^(cos0i) A *(cos0 2 )' lj , (17) 



where cos 9i is the cosine of f*j with the z-axis. The de- 
termination of the expansion coefficients gj is straightfor- 
ward but a little cumbersome. Except for an insignificant 
overall factor, the gj are rational numbers with not too 
large denominators, which allowed us to compute them 
numerically and ensure that they were accurate to all 
digits of arithmetic precision. 

For the angular integration we change variables to 
cos 0i, ipi, cos #12 and tpi2, where tp% is the azimuthal an- 
gle of fi, cos6>i2 = f\ ■ ?2, and <^i2 is the azimuthal angle 
of r*2 with respect to the axis fi. Integrations can then 
be performed over all angular variables except for cos 6*12, 
which is expressed as cos 6*12 = (rf + r\ — r\^) / (2r\r2) ■ 
The remaining three-dimensional integrals have the gen- 
eral form 



ar 1 -f3r 2 



(18) 



The integrals are non-trivial, because r\, and r\% are 
connected by the triangular inequality. We compute 
them by a strictly positive and therefore numerically sta- 
ble recurrence formula |25| . 

Loss of accuracy is due to the expansions of Gu and 
cos #12 , which at high L become very lengthy with up 
to several hundred terms. Beyond L sa 5 we needed to 
use Fortran REAL* 16 arithmetic, which provided stable 
results up to L w 12. 

Even when no accuracy is lost in the matrix elements 
the metrical matrix S has in general a very poor condition 
number which quickly exceeds regular machine precision 
14 decimal digits) or even REAL* 16 precision (« 32 
decimal digits). This may lead to uncontrolled errors. 
We remove the problem by first rescaling the basis such 
that the diagonal elements of the metrical matrix are 
Su = 1, after which the maximum size of the eigenvalues 
of S was limited to < 100. We diagonalize S and remove 
eigenvectors with eigenvalues smaller than a threshold e. 
We found e = 10~ n to be suitable for Fortran REAL*8 
arithmetic. In REAL* 16 we could use e = 10~ 30 , but 
results were found to depend very weakly on the thresh- 
old. The remaining eigenvectors are normalized with 
respect to S 



After transforming all matrices to the orthonormal basis 
{!&)} precision can be lowered to REAL*8 to save on 
storage and computation time. 



E. Time propagation 

For the time propagation we make one more transfor- 
mation to the "atomic basis" that diagonalizes the 
complex scaled atomic Hamiltonian 



(r)i\Ho,e\rij) = E i8 



(20) 



where 



-(Ai + A 2 ) 



Ze 



Ze 



ri2 



(21) 



In dipole approximation the laser field only connects an- 
gular momenta differing by \L — L'\ = 1. The resulting 
structure of the overall Hamiltonian is depicted in Fig. ||. 

At the laser intensities and frequencies in our calcula- 
tions the diagonal terms still make the dominant contri- 
bution to the time evolution. Therefore we time integrate 
in the "interaction picture" 



z{t) = a(t) + b(t) 



*37 a W = H ,ea(t) 
at 



H e (t)-H ,e a(t) + H e (t)b(t). 



(22) 
(23) 

(24) 



Here c(t) denotes the coefficient vector and 77 ,e and He 
are the operator matrices with respect to the {|?7i)} basis. 
For the decomposition of c(t) into a(t) + b(t) we set at 
the beginning of each time step 

a(t ) = c(to), b(t ) = 0. 

Since H g is diagonal in our basis, the solution of Eq. ( |2^ ) 
is trivial. Eq. (24) is solved by a seven stage 6th order 
explicit Runge-Kutta method (Butcher's method, given 
in Ref. [^6|), which we found to be more CPU time effi- 
cient than lower order methods. There were no problems 
with numerical stability, as we verified by comparing with 
lower order methods at a few parameter points. The time 
step was automatically adapted by comparing every two 
integration steps with a single double step size integra- 
tion. The typical number of time steps was about 400 
per optical cycle, which for our 7 stage method means 
about 3000 matrix-vector multiplies per cycle. Compu- 
tation times for the shorter pulses were about 1 hour on 
a 500 MHz DEC/ Alpha work station. 



(19) 
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III. RESULTS 
A. Bound and doubly excited states 

Before we discuss the results of the time propagation, 
we want to list the energies and widths of the most im- 
portant bound and doubly excited states of He and H~ 
as obtained with the above basis. Table || compares the 
values of the first few bound state energies of He that 
we obtained with the basis used in the time propaga- 
tion with reference values from literature. Several of our 
values are lower than the variational upper bounds from 
literature, but this does not indicate greater accuracy, 
since due to complex scaling our values are not upper 
bounds. Our accuracy is ~ 10~ 8 a.it. for most energies 
given. The basis sizes are about 300 for each angular 
momentum. The functions are counted after removal of 
near-singular vectors from the basis (cf. section |5 D| ), 
which is the number relevant for the time propagation. 
As an example, we listed in Table | above the basis set 
for L = 2. It has been shown that literature values can 
be exactly reproduced with bases of size ~ 150 that are 
optimized for each state fL9| , pO| ]. 

Table III gives lowest few doubly excited states of He 
from the time-propagation. The states are labelled by 
the radial quantum numbers of the dominant single elec- 
tron contributions n\ and n-i- Accuracies are of the order 
10~ 4 a.u. for the energies and widths. At the very small 
widths the relative accuracies can become poor. 

Finally Table IV summarizes the bound and doubly 
excited state energies of H~. We present the energies 
as they appear in the time propagation as well as values 
obtained in larger bases with state specifically adapted 
exponents a s and [3 S . Judging from the convergence be- 
havior, we believe that our state-specific values are accu- 
rate to all except possibly the last digit quoted. One sees 
that some of them are more accurate than the literature 
values quoted. A special case is the shape resonance in 
L = 1 just above the H (n = 2) threshold. It requires 
a minimum scaling angle of 9 f=a 0.25 to become mani- 
fest as an isolated eigenvalue of the scaled Hamiltonian. 
With the smaller angle used in the time propagation the 
resonance state cannot be distinguished from the approx- 
imate continuous states that surround it. 



B. Excitation and ionization of He 

Fig. U shows the probability of excitation and ioniza- 
tion of He by cos 2 -shaped pulses of duration T — 157 a.u. 
and peak intensity I = 0.00423 a.u. = 2.97 x 10 14 W/ cm 2 . 
Note that in Ref. Q the conversion to SI units was 
too small by a factor of 2. The frequencies cover the 
range from just above the five photon ionization thresh- 
old to well above the single ionization threshold. Below 
each of the thresholds one clearly distinguishes the en- 
hancement of bound state excitation due to resonances. 



The peaks below the one photon and two photon thresh- 
olds are due to resonances with the lowest P-state (en- 
ergy = —2.12384) and the lowest excited 5-state (en- 
ergy = —2.14597), respectively. The resonances below 
the lower thresholds are not well separated due to the 
spectral width of the pulse. The pronounced dip in the 
bound state excitation at w = 0.72 is due to a Rabi- 
like oscillation. When one increases the pulse duration, 
the minimum disappears completely and reappears at a 
pulse duration of ~ 270. A third minimum appears at 
T = 380 (Fig. |). The oscillation period roughly cor- 
responds to the Rabi period at the given parameters. 
The effect strongly depends on the pulse shape. Fig. |] 
compares excitation by cos 2 -pulses with excitation by 
Gaussian pulses. The pulse energy is 0.5 a.u. for both 
pulse shapes and the widths of the vector potential en- 
velopes was T — 157 for the cos 2 -pulse and T = 92 for 
the Gaussian pulse. With this choice the envelopes have 
the same width at 1/e of the maximum for both pulse 
shapes. One can see that the Rabi-like oscillations occur 
for both shapes, but they are more pronounced with the 
cos 2 -pulses. 

One can also distinguish resonant enhancement of ion- 
ization, manifested by the coincidence of bound-state ex- 
citation peaks with ionization peaks. Most pronounced 
are the peaks at the three-photon resonance with the 
lowest P-state at frequency u> — 0.26 and the two-photon 
resonance with the lowest excited 5-state at u> = 0.38; 
another resonance with the lowest P-state is hidden in 
the slope to the two-photon ionization threshold. In 
spite of the massive bound state excitation at the single- 
photon resonance at uj = 0.78, weak coupling of the P- 
state to the continuum leads only to a slight bump in the 
ionization rate at that frequency. 

The lowest frequency of 0.2 in the figures is somewhat 
above the the popular frequency of 0.1837 of a 248nm 
wavelength laser. In Table [v] we give the ionization 
yield at wave length 248 nm by a cos 2 pulse of dura- 
tion T = 40 cycles = 32.8 fs, which has a half width 
of 20 optical cycles. The pulse parameters coincide with 
the ones used in Ref. |3lf | for the comparison with experi- 
ment. Our results are systematically lower by about 50% 
compared to Ref. |$(|. We estimate that our results are 
converged to 10% accuracy or better (see below). Consid- 
ering that we cover intensities all the way into the onset 
of saturation, the agreement is nevertheless satisfactory. 

Recently another calculation for the same pulse has 
been published |37[ • Those values sizably oscillate around 
ours and also around the results of Ref. |56|. The dis- 
agreement is particularly evident at the lower intensities, 
where convergence problems should be smaller. It ap- 
pears therefore that calculations of Ref. |37j remain rel- 
atively far from convergence. A possible reason for the 
lack of convergence is the absence of doubly excited states 
from the basis of Ref. p7| ], although the velocity gauge 
used in that calculation always introduces at least virtual 
doubly excited states (cf. section II B). 

The experimental number given in Ref. pq| is at least 
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one order of magnitude smaller than all theoretical re- 
sults and none of the calculations falls within the quite 
large upper error margin of the experiment. 

Ionization is predominantly a single electron effect and 
all observations made above are qualitatively the same in 
a single electron atom. Double excitation in turn is a gen- 
uine two electron effect. Fig. |^ shows the population of 
the lowest doubly excited states with angular momentum 
L = 0, 1, 2, and 3 after passage of a cos 2 -shaped pulse 
of duration T = 157. The narrower peaks are due to 
multiphoton resonances between the ground and the re- 
spective doubly excited states with 2 to 6 photons. The 
population of the states with even L are enhanced in the 
range u> = 0.6 ~ 0.9. This enhancement is almost exactly 
proportional to the 1 P° bound state population and the 
ratio between the populations was found to be propor- 
tional to the laser intensity. We therefore interpret it as 
a far off-resonant single photon transition from the P 
to the S and D states brought about by the large band 
width of the pulse. As with the Rabi oscillation of sin- 
gle excitation, the pulse shape is quite important, since 
with Gaussian pulses the phenomenon nearly disappears. 
Similarly, longer pulse duration suppresses the effect. 



C. Excitation and electron detachment from H 

The binding energy of the only bound state of H~ is 
0.02775 a. u.. Therefore much lower laser intensities and 
frequencies lead to total electron detachment from H~ . 

Fig. U shows the photodetachment at a laser frequency 
of to = 0.03 for intensities between 1 x lO^W/cm 2 and 
8 x 10 n W/cm 2 . At about 2 x 10 n W/cm 2 the single 
photon ionization threshold rises above ui = 0.03 due to 
the AC Stark shift and roughly in the same intensity re- 
gion two photon ionization becomes dominant |5| . In our 
time-dependent calculations we can distinguish a bend 
in the detachment yield, which moves towards expected 
intensity of 2 x 10 11 W/cm 2 with increasing pulse dura- 
tion. The dashed lines are obtained by integrating the 
intensity dependent detachment rates r(J) from Ref. || 
with the pulse shapes used in our calculations. Satura- 
tion effects in the final detachment yield Y(t — oo) are 
included according to the equation 



Y(t) 



dt'T[I(t')} [l-Y(t% 



As the instantaneous intensity we defined 

J(t) := [E Q h(t)] 2 /2, E Q = coA . (26) 

This definition neglects terms with the time-derivative of 
the envelope, whose contributions do not visibly change 
the comparison. The agreement for the two longer pulses 
is quite good, except that the bend in the yield is some- 
what more pronounced for the integrated rates. With 
the shorter pulses the integrated rate underestimates the 



true detachment. Assuming the applicability of the rate 
concept also at the short pulses, the spectral width of the 
pulse could be sufficient to explain the higher yield: with 
a broad pulse the closure of the single photon ionization 
channel is moved to higher intensities, thus effectively en- 
hancing the yield. For a quantification of this hypothesis 
rates T(I) for frequencies other than uj = 0.03 would be 
required. 

For the double excitation of H~ one needs frequen- 
cies in the visible to low UV. Fig. ^ shows the excitation 
of the relatively long-lived lowest autoionizing P state 
around the resonant frequency of u) = 0.134 and intensity 
10 13 VF/ cm 2 . The resonance is shifted from its field-free 
position by about 0.001 for the 10 cycle pulses. Longer 
pulses of 20 cycles have a larger shift of ~ 0.002. In the 
limit of a constant laser field an AC-Stark shift of the 
ground state by ~ 0.005 is expected by perturbatively 
extrapolating the data from Ref. |5j 



D. Harmonic generation 

By Eq. (||) we calculate the expectation value of the 
dipole d(t) as a function of time, from which we obtain 
the acceleration of the dipole by numerical differentia- 
tion. It would be more desirable to directly calculate the 

expectation values of the acceleration of the dipole d(t), 
but unfortunately in the matrix elements of that operator 
integrals of the form ( |l8|) with negative powers of r\ or 
ri arise, which can only be calculated with considerable 
extra numerical effort (cf. p5|). 

The expectation values d(t) are sensitive to the wave 
function at larger distances. To obtain accurate results 
we therefore added extra basis functions with smaller ex- 
ponents to cover a longer range in r\ and r%- By this 
enlargement of the basis we were able to obtain satis- 
factory accuracies of < 10% up to the 5th harmonic for 
frequencies in the range 0.34 to 0.42 a. u. and intensity 
I = 2 x 10 14 W/cm 2 . The frequencies include a strong 
resonance with the lowest excited S-state that greatly en- 
hances harmonic generation and leads to the dominance 
of the 3rd harmonic over the 1st. Figure || summarizes 
the results obtained with a pulse duration of 40 optical 
cycles. Harmonics up to order 7 can be distinguished and 
the resonant enhancement at u> = 0.38 is manifest. For 



(25) reference we include Table VI with the peak heights 



E. Accuracy of the results 

The limited expansion length for the wave function in- 
troduces the dominant error into our calculations. Now 
we study the effect of the truncations with respect to 
total angular momentum and radial basis functions. 

Fig. |^ compares the results of calculations with L max = 
5, 6, and 7 at peak intensity 2.97 x 10 14 IV/ cm 2 and fre- 
quencies between 0.2 and 0.4. Only at the lowest frequen- 
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cies there is a distinguishable effect on ionization with 
a relative difference between the calculations of < 1%- 
Double excitation is more sensitive to angular momentum 
truncation, since it is a much higher order process, but 
still the relative error rises to only about 10%. Assum- 
ing exponential convergence, we can conclude that with 
L max also for double excitation the error due to angular 
momentum truncation is < 1% at the given parameters. 

Convergence of the expansion in the internal coordi- 
nates n, r2, r\2 is more difficult to investigate, since the 
basis rapidly grows when one increases the admissible 
powers of ri, r2, t\i- Complex scaling provides an indi- 
rect estimate of the accuracy of the internal expansion. 
For an infinite basis, the results are independent of 9. 
Any dependence on 9 must therefore be ascribed to the 
basis truncation. In practice, there is only a limited range 
of 9 where the results vary little with 8. When 9 is too 
small, the outgoing waves are only weakly damped and 
one needs to describe a long oscillatory tail in the wave 
function, for which an L 2 -expansion converges slowly. 
When on the other hand 9 is too large the complex scaled 
bound state wave functions have increasingly oscillatory 
character, which again is not well reproduced by the fi- 
nite basis (cf . appendix) . Fig. [l0| shows ionization and 
double excitation obtained with different scaling angles 
at laser frequency 0.39 a.u.. One distinguishes a range of 
9 where the results are quite stable. The variation inside 
this range is of the same size as the variation when we 
increase the number of basis functions by a factor ~ 2, 
which supports our use of the variation with 9 for an 
accuracy estimate. For the given parameters ionization 
varies less than 0.2%. The accuracies for double excita- 
tion are only slightly lower. 

At frequency uj = 0.3 the variation of ionization still re- 
mains within the 1% range, but double excitation varies 
by about 10%, which indicates that we approach the lim- 
its of numerically reliable results for multiphoton double 
excitation. 



IV. SUMMARY AND CONCLUSIONS 

The method introduced in this paper allows the nu- 
merical integration of the complete Schrodinger equation 
of a two-electron atom in a strong laser pulse with re- 
alistic frequency, intensity, and duration. This has been 
demonstrated on the systems of Helium and H ~ with fre- 
quencies ranging from infrared to ultraviolet, intensities 
up to 10 15 W/cm 2 and pulses as long as 160 fs. 

The first important constituent of our method is an 
explicitly correlated basis set expansion, which allows to 
reproduce the atomic structure to essentially any desired 
accuracy at moderate expansion length. This includes 
positions and widths of doubly excited states, where we 
obtain accuracies that rival and in some cases exceed lit- 
erature values. 

The second ingredient of our method, complex scaling, 



was formally introduced as a method of imposing strictly 
outgoing boundary conditions. We cannot at present give 
a mathematically rigorous theory for the use of complex 
scaling in a time-dependent calculation, but we provide 
heuristic arguments and numerical evidence that it in- 
deed is equivalent to the regular Schrodinger equation 
with strictly outgoing boundary conditions. The techni- 
cal advantage of complex scaling is that the expansion 
length remains short. 

Convergence was investigated for the whole range of 
parameters and indicates accuracies between fractions of 
a percent at higher laser frequencies and at least 10% for 
the majority of the data. Only at the seventh harmonic 
peak in our example and for double excitation at lowest 
laser frequencies accuracies remain unsatisfactory. 

The numerically most challenging combination of pa- 
rameters was used for Helium exposed to pulses of du- 
ration ~ 32 fs at the wave length of 248 nm and peak 
intensities up to 10 15 W/cm 2 , where we reach an accu- 
racy of the ionization yield of about 10%. A previous 
calculation for the same pulses (36| qualitatively agrees 
with ours but exceeds our result by about 50% at the 
highest intensity. A more recent calculation j3^j] deviates 
more strongly also at the lower intensities. 

Another comparison with existing theoretical work 
could be performed with a Floquet calculation for the 
electron detachment from H~ . Quite satisfactory agree- 
ment was found for pulses of at least 16 optical cycles. 
Shorter pulses may not be expected to compare well with 
a Floquet calculation for constant intensity. 

We believe that our results should serve as a bench- 
mark for future calculations of two-electron systems. 
This refers to both, complete two-electron calculations 
as well as model calculations. It may be expected that 
for harmonic generation and ionization a single-electron 
description will be found to be satisfactory in a range of 
parameters. 

If necessary, our method allows extensions in several 
directions. To extend the range of accessible parameters 
is predominantly a question of more computer power, al- 
though also technical modifications in the calculation of 
the matrix elements are required to control the loss of 
accuracy. A second obvious extension is the introduc- 
tion of an electronic core to model effective two-electron 
atoms like Mg. Finally, the interpretation of the com- 
plex scaled wave function adopted in this paper, which 
takes the back-scaled wave function as an approxima- 
tion to the regular solution with outgoing wave boundary 
conditions, suggests that electron spectra can be deter- 
mined. Whether this is numerically feasible and practical 
remains to be investigated. 
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APPENDIX A: 

Here we deduce the relation between outgoing wave 
boundary conditions and complex scaling. We first as- 
sume that one can write the radial wave function in the 
form 



¥(r;t) 



dkc(k,t)$ k (r) 



(Al) 



where the have the asymptotic behavior ^ < 
^ solves the Schrodinger equation 



i-V(r;t) = H(r;t)*(r;t). 



■ikr 



and 



(A2) 



For the sake of brevity we have omitted the part of the 
expansion with square-integrable functions. (For the 
Coulomb potential the asymptotic behavior is more pre- 
cisely ~ exp(i/cr— i lrx2| fc[r/ fc).) In terms of the expansion 
coefficients c(k;t) Eq. (A2) can be written as 



i Jt c(k;t) = 



dk'h(k,k';t)c(k';t). 



(A3) 



Outgoing boundary conditions mean that one solves 
Eq. (IX|) restricted to k > 0: 



dk'h(k,k';t)c+(k';t). 



(A4) 



The time dependent wave function with outgoing bound- 
ary conditions is then 



dk'c+(k';t)<f> k ,(r) 



(A5) 



In order to relate ^f(r,t) to the complex scaled wave 
function we must assume that ty(r,t) and its time- 
derivative d^S (r, t) I dt are analytic functions of r for any 
analytic initial state ^S(r,t — 0). This assumption is non- 
trivial: for example, it is known to be valid for the Hamil- 
tonian of the field-free hydrogen atom, while it is obvi- 
ously violated for potentials that are non-differentiable 
at any point other than r — 0. Under this assumption the 
analytically continued wave function ty(j]r; t), Im(rj) > 
solves the "complex scaled" Schrodinger equation 



/oc 
dkc v (k,t)$ k (rir), 
-OO 

do not depend on i] 



; (A,t) = c(*;,t). 



(A7) 



(A8) 



From that it follows that also the kernel of the time- 
integration h(k, k'; t) does not depend on -q. We see that 
the two equations ( |Ag ) and (A6) describe exactly the 
same dynamics. The important difference is that due 
to the asymptotic behavio r of the expansion functions 
&k{w) ~ exp(z/c?/r) in Eq. (A6) we can distinguish k > 
from k < by the norm: "ingoing waves" k < grow 
exponentially, while "outgoing waves" k > become 
square-integrable. Because of the exponential divergence 
any function 



dka(k)^ k (rjr) 



(A9) 



will diverge, if /_ dk\a(k)\ 2 > 0. Consequently, if 
^{•qr;t) contains ingoing waves, it will not be square in- 
tegrable. To obtain a solution with ou tgoing waves only, 
we solve the differential equation ( |A6| ) restricted to the 
space of square integrable functions | \^ + (r]r; t)\ \ < oc. 
Since h(k,k';t) does not depend on 77 the coefficients of 
the expansion 



dk'c + (k';t)<f> k ,(rir) 



(A10) 



are the same as in Eq. ( |A5| ) and the outgoing wave solu- 
tion is obtained by substituting r\r w ith r. 

As initial condition for Eq. (A6) we use a field free 
bound state. For the class of "dilation analytic" po- 
tentials, which include the Coulomb potential, the bound 
state functions are known to be analytic functions of r. 
This is trivial to verify for the complex scaled Hamilto- 
nian of the two-body Coulomb problem 

-— ^2 A - — ) $i(r)r) = Ei$i(rfr), (All) 
ay ijr J 

where 3>i(r) is a hydrogenic bound state function. For 
computation it is useful to keep in mind that, for exam- 
ple, the radial ground state eigenfunction r;rexp(— rjr) 
becomes increasingly oscillatory with increasing 7m(?y), 
which is opposite to the outgoing waves, where larger 
Im(ji) causes stronger damping. Oscillatory functions 
are generally more difficult to represent numerically and 
the choice of 77 = e %e will depend on whether the outgoing 
wave or the bound state part of ^ v are more important. 



H{r,r;t)^(r]r;t), 



(A6) 



which is the analogue of equation (||) for a single ra- 
dial coordinate. If we further assume that the expansion 
functions <fr k are analytic, the expansion coefficients c v 
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TABLE I. Basis set for L = 2 used in time propaga- 
tion. "Size" denotes the number of basis functions. For 
the time-propagation near-linearly dependent vectors are re- 
moved. 
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6 6 
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6 


139 


total 










416 


total in time-propagation 
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TABLE II. 


Bound state energies of He: (a) 


calculated 


with the basis used for time propagation, (b) literature values. 


(a) 




(b) 
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L = 5 




L = 7 










-2.0138890346 




-2.0078125124 










-2.0102041827 




-2.0061728489 










-2.0078125737 
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TABLE III. Doubly excited states of He for L = through 
3. (a) values obtained with the basis used in the time propa- 
gation, (b) literature values from references [29] (L < 2) and 
[28] (L > 2). 

(ni,n 2 ) (a) (b) 
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TABLE IV. Bound and doubly excited states of H~ for TABLE V. Ionization yield for laser wave length 248 nm 

L = through 3. (a) values obtained with the basis used in and pulse duration T = 40 optical cycles as a function of peak 

the time propagation, (b) values with a basis optimized for intensity I. The literature values are obtained by converting 

each state, and literature values. The values (b) are estimated the generalized cross sections from Table 2 in Ref. [36] with 

to be converged to all digits given except for the last. the help of Eq. (16) in that reference. 



Energy 

L = 

-0.52775101689 
-0.5277510165443 





Width 




1.7324xl0~ 3 
1.7332 xl0~ 3 
1.731 xl0~ 3 

12 xlO" 5 
9.02 xlO" 5 
8.985 xl0~ 5 


Reference 

present, (a) 
@ 


I(W/cm 2 ) 
2 x 10 14 
2.5 x 10 14 
5 x 10 14 
1 x 10 15 


present Ref. [|36[ 
7.06 x 10 -3 8.13 x 10~ 3 
0.00105 0.00148 
0.043 0.069 
0.18 0.33 


-U.14877D4 
-0.1487762 
-0.1487765 

-0.12605 

-0.1260199 

-0.12601965 






present, 
present, 


present, 
present, 

n 


(a) 
(b) 

(a) 
(b) 


TABLE VI. Relative peak heights of harmonic generation 
by cos 2 -pulses of T = 40 optical cycles and peak intensity 
2 x W 14 W/cm 2 for three different fundamental frequencies to 
The accuracies are 10% up to the 5th harmonic and of the 
order 50% for the 7th harmonic. 


-0.069006 
-0.069006 






1 41 99 x 10 -3 
1.4184 xlO" 3 


present, 
present, 


(a) 
(b) 


LO 1 


Harmonic order 
3 5 7 


-0.05615 
-0.0561434 






23 xlO" 5 
8.8 xlO -5 


present, 
present, 


(a) 
(b) 


0.34 0.85 
0.38 1 
0.42 1.46 


0.25 4.8 x 10" b 4.4 x 10" s 
3.9 3.2 x 10~ 4 3.4 x 10~ 7 
0.11 2.1 x 10" 6 


-0.12604986 
-0.12604986 
-0.1260495 


L = 


1 


1.36 x 10" 6 
1 36 x in -6 
1.165 x 10" 6 


present, 
present, 

@ 


(a) 
(b) 






*) 

-0.1243856 
-0.12436 






l .U A 1U 

6.9 x 10" 4 


present, 
present, 

@ 


(a) 
(b) 






-0.062708 
-0.062716 
-0.06871675 






1.17 x 10" 3 
1.19 x 10" 3 
1 1 Q1 4 x 1 -3 


present, 
present, 

@ 


(a) 
(b) 






-0.058586 

-0.0585718 

-0.05857181 






< 10" 5 
8.988 x 10~ B 
8 Q86 x 1 -6 


present, 
present, 


(a) 
(b) 






-0.127937 
-0.127937 
-0.12794175 


L = 


2 


3 19 x10~ 4 
3.12 xlO" 4 
3.1625xl0 -4 


present, 
present, 


(a) 
(b) 






-0.065954 

-0.0659531 

-0.0659533 






1.D04 X 1U 

1.6576 xlO" 3 
1.6581 xl0~ 3 


present, 
present, 

@ 


(a) 
(b) 






-0.056834 
-0.0568294 






2.8 xlO -4 
2.5302 xlO -4 


present, 
present, 


(a) 
(b) 






-0.056564 
-0.05655875 


L = 


3 


3.54 xl0~ 3 
5.00 xlO -3 


present, 
present, 


(a) 
(b) 






*) state cannot be distinguished from surrounding 


continuum 


at 







the complex scaling angle of 8 = 0.22 used in time propagation 
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FIG. 1. The Hamiltonian matrix in the atomic basis 
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CD 
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0.4 0.6 0.8 1.0 
Frequency (a.u.) 
FIG. 2. Excitation and ionization of He by a cos 2 -shaped 
pulse of duration 3.8 fs and peak intensity 2.97 x W 14 W/crn 2 
as a function of frequency. Solid line: ionization, dashed line: 
ionization plus bound state excitation. The arrows labelled 
by n = 1, 2, 3, and 4 indicate n-photon ionization thresholds. 
The dip at frequency 0.72 is due to a Rabi oscillation. 
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FIG. 3. Bound state excitation as a function of pulse 
duration at frequency u> = 0.72. Solid line: total population of 
excited bound states. Dashed line: ionization. The distance 
between the minima is roughly the Rabi period. 
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FIG. 4. Bound state excitation for Gaussian and cos 2 
pulses as a function of frequency. The pulse energy is 0.5a.u. 
for both shapes, the pulse widths are T = 157 for cos 2 pulses 
(solid line) and T = 92 for Gaussian pulses (dashed line) . 
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FIG. 5. Populations of the lowest doubly excited states 
with L = 0,1, and 2 after the passage of a pulse with 
T = 3.8 fs with peak intensity 2.97 x 10 14 W/cm 2 . Dashed 
lines: population of the lowest bound state with symmetry 
1 P° x!0~ 4 . 
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FIG. 6. Ionization of H~ by cos 2 -shaped pulses with laser 
frequency w = 0.03 as a function of intensity for pulse dura- 
tions T = 8, 16, 24, and 32 optical cycles. The bend in the 
ionization yield is due to the closure of single-photon ion- 
ization channel by the AC-Stark shift and the transition to 
two-photon ionization. Dashed lines: time-integrated Flo- 
quet rates from Ref. [5] 




0.125 0.130 0.135 0A40 

Frequency (a.u.) 
FIG. 7. Population of the lowest doubly excited P-state 
of H~ by cos 2 -shaped pulses as a function of frequency. Peak 
intensity is 10 13 W/cm 2 and pulse duration duration T = 10 
optical cycles (solid line) and T — 20 (dashed line). 
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FIG. 8. Harmonic generation on Helium with a cos 2 -pulse 
of duration T = 40 optical cycles and peak intensity 
2 x 10 14 W/cm 2 . The fundamental frequencies are 0.34 
(dot-dashed line), 0.38 (solid line) and 0.42 a.u. (dashed 
line), respectively. The frequency 0.38 is two-photon reso- 
nant with the lowest excited S-state. The structure around 
ui — 0.8 originates from bound state excitations. 
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FIG. 9. Convergence of ionization and total double exci- 
tation with the maximum angular momentum L max . Dashed 
line: relative difference between calculations with Z/ max = 5 
and 7; solid line: relative difference between L max = 6 and 7. 
Intensity = 2.97 x W 14 W/cm 2 , pulse duration T = 157. 
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FIG. 10. Dependence of ionization and double excitation 
on the complex scaling angle 9. The dotted lines indicate esti- 



mated errors. Intensity = 2.97 x 10 
and pulse duration T = 157. 



W/cm , frequency=0.4 
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